//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// cl_fit_mv - Ajusta curva de luz de binarias a dados observados utilizando
//			   Simplex. 
// 
// Definicao das funcoes da classe PReader que le os parametros de entrada
// para programa de ajuste da curva de luz. Essa classe le e armazena os parametros
// de uma forma amigavel ao programa de modelagem de curva de luz. Alem disso, 
// a classe define a interface para passar os dados aos programas de modelagem
// atravez da funcao orgpar().
// 
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

//#include "cl_fit_input_funcs.hpp"
//#include "cl_fit_data_funcs.hpp"
//#include <exception>
#include "cl_fit_mc_v0_1.h"
#include "dataFileReader_v0.1.1/dataFileReader.h"
#include "light_curve.h"

double chi_square(const gsl_vector *par, void * GSLSimplex);

using namespace std;

//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Definindo construtor

GSL_Simplex::GSL_Simplex(char* arg1, char* arg2):PReader(arg1),DataReader(arg2){
	//
	//--------------------------------------------------------
	// Definindo parametros a ser ajustados. Defino uma forma padrao.
	//
	_FitPar = new double * [ifdnp()+nspot()]; // Parametros a serem ajustados
	_PrivParams = new double * [ifdnp()+nspot()]; // Quais parametros sao privados (diferentes para cada curva de luz)
	
	//
	//--------------------------------------------------------
	// Aproveito heranca para ler propriedades dos arquivos de entrada :)
	//
	for(int i = 0; i < ifdnp(); i++){
		
		_FitPar[i] = new double [ifdp(i)];
		_PrivParams[i] = new double [ifdp(i)];
		
		for(int j = 0; j < ifdp(i); j++){
			_FitPar[i][j] = 0.0; // read.readColumn<double>(j+1);
			_PrivParams[i][j] = 0.0;
		}
		
	}
	
	for(int i = ifdnp(); i < ifdnp()+nspot(); i++){
		
		_FitPar[i] = new double [ifdsp()];
		_PrivParams[i] = new double [ifdsp()];
		for(int j = 0; j < ifdsp(); j++){
			_FitPar[i][j] = 0.0; 
			_PrivParams[i][j] = 0.0;
		}
	}
	
	const int NDefFitPar = 6; // Numero de parametros ajustados default
	const int NDefPrivPar = 2; // Numero de parametros privados default (f1 e f2)
	
	_FitPar_index = new int* [2];
	_PrivPar_index = new int* [2];
	
	_FitPar_index[0] = new int [NDefFitPar];
	_FitPar_index[1] = new int [NDefFitPar];
	
	_PrivPar_index[0] = new int [NDefPrivPar];
	_PrivPar_index[1] = new int [NDefPrivPar];
	
	// inc, M2, f2, f_fac2, f1, f_fac1
	int i_DefFitPar[] = {1,1,2,2,3,3}; // linha do vetor de parametros a ser ajustado
	int j_DefFitPar[] = {0,2,0,5,0,5}; //	coluna do vetor de parametros a ser ajustado
	
	int i_DefPrivPar[] = {2,3}; // linha do vetor de parametros a ser ajustado
	int j_DefPrivPar[] = {0,0}; //	coluna do vetor de parametros a ser ajustado
	
	//
	//--------------------------------------------------------
	// Seto os parametros a serem ajustados e guardo os indices
	//
	for (int i = 0; i < NDefFitPar; ++i)
	{
		_FitPar[i_DefFitPar[i]][j_DefFitPar[i]] = 1;
		_FitPar_index[0][i] = i_DefFitPar[i];
		_FitPar_index[1][i] = j_DefFitPar[i];
		
		if( i < NDefPrivPar ){
			_PrivParams[i_DefPrivPar[i]][j_DefPrivPar[i]] = 1;
			_PrivPar_index[0][i] = i_DefPrivPar[i];
			_PrivPar_index[1][i] = j_DefPrivPar[i];
		}
	}

	_NFitPar = NDefFitPar; //guardo numero de parametros a ajustar
	_NPrivPar = NDefPrivPar; //guardo numero de parametros privados a cada curva de luza
	
	//
	//--------------------------------------------------------
	// inicializando parametros da reconstrucao
	//
	
	_nIter = 100; // numero de iteracoes antes de finalizar programa. As interacoes sao contadas a cada troca de temperatura.
	_iter = 0;
	
	_defautSSize = 0.25;
	
	_StepSize = new double [_NFitPar]; // Tamanho do passo eh 10 por cento do valor do parametro
	
	for (int i = 0; i < NDefFitPar; ++i)
	{
		if( fabs(getpar(i_DefFitPar[i], j_DefFitPar[i]) ) > 0.0 )  {
			_StepSize[i] = fabs(getpar(i_DefFitPar[i], j_DefFitPar[i]) ) * _defautSSize;
		}
		else {
			_StepSize[i] = _defautSSize; 
		}
	}
	
	//
	//--------------------------------------------------------
	// inicializando funcoes da GSL
	//
	
	T =	gsl_multimin_fminimizer_nmsimplex2;
	
	x = gsl_vector_alloc(_NFitPar+_NPrivPar*(nData()-1));
	ss = gsl_vector_alloc(_NFitPar+_NPrivPar*(nData()-1));

	for (int i = 0, j = 0; i < _NFitPar; ++i){

		if (_PrivParams[_FitPar_index[0][i]][_FitPar_index[1][i]] > 0.0){
			// Parametro eh privado de cada curva de luz
			for (int k = 0; k < nData(); ++k) {
				gsl_vector_set(x, i+k, getpar(i_DefFitPar[i], j_DefFitPar[i]));
				gsl_vector_set(ss,i+k, _StepSize[i]); // tamanho do passo p/ funcao de minimizacao
			
			}
			j += nData()-1;
			
		}
		else {
			gsl_vector_set(x, i+j, getpar(i_DefFitPar[i], j_DefFitPar[i])); // Parametro para funcao de minimizacao
			gsl_vector_set(ss, i+j, _StepSize[i]); // tamanho do passo p/ funcao de minimizacao
		}
	

	}
	
	minex_func.n = _NFitPar+_NPrivPar*(nData()-1);
	minex_func.f = &chi_square;
	minex_func.params = (void *)this;
		
	// Raiz do nome do arquivo de saida
	_OutputRootFilename = "simansolve";

}

// Fim do contrutor
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Fucao que mostra parametros a serem ajustados

void inline GSL_Simplex::show_FitPar(ostream& out)
{
	cout << "# - Fitting parameters:\n";
	
	for(int i = 0; i < ifdnp(); i++)
	{
		out << "# ";
		
		for (int j = 0; j < ifdp(i); j++) if(_FitPar[i][j]) out << pardef(i,j) << getpar(i,j) << " | ";
		
		out << endl;
	}
	
	for(int i = ifdnp(); i < ifdnp()+nspot(); i++){
		
		out << "# ";
		
		for(int j = 0; j < ifdsp(); j++) if(_FitPar[i][j]) out << pardef(ifdnp(),j) << getpar(i,j)  << " | ";
		
		out << endl;
		
	}
	
	
}

// Fim fucao que mostra parametros a serem ajustados
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Fucao que muda parametros a serem ajustados

void inline GSL_Simplex::FitPar(const int npar, int* i_pos,int* j_pos, const int npriv, int* i_priv,int* j_priv)
{
	// Zerar escolha de parametros
	
	for(int i = 0; i < ifdnp(); i++){
		for (int j = 0; j < ifdp(i); j++) {
			_FitPar[i][j] = 0;
			_PrivParams[i][j] = 0;
		}

	}
	
	
	for(int i = ifdnp(); i < ifdnp()+nspot(); i++){
		for(int j = 0; j < ifdsp(); j++){
			_FitPar[i][j]  = 0;
			_PrivParams[i][j] = 0;
		}

	}
	
	
	// Ajustando tamanho do passo
	delete _StepSize;
	_StepSize = new double [npar]; 
	
	delete _FitPar_index[0];
	delete _FitPar_index[1];
	
	_FitPar_index[0] = new int [npar];
	_FitPar_index[1] = new int [npar];
	
	delete _PrivPar_index[0];
	delete _PrivPar_index[1];	
	
	_PrivPar_index[0] = new int [npriv];
	_PrivPar_index[1] = new int [npriv];
	
	// Recelecionando parametros
	for (int i = 0; i < npar; ++i) {
		_FitPar[*i_pos][*j_pos] = 1;
		_FitPar_index[0][i] = *i_pos;
		_FitPar_index[1][i] = *j_pos;
		
		if( fabs(getpar(*i_pos, *j_pos) ) > 0.0 )  
			_StepSize[i] = fabs(getpar(*i_pos, *j_pos) ) * _defautSSize;
		else _StepSize[i] = _defautSSize; 
		
		if(i < npriv){
			_PrivParams[*i_pos][*j_pos] = 1;
			_PrivPar_index[0][i] = *i_priv;
			_PrivPar_index[1][i] = *j_priv;
			++i_priv;
			++j_priv;
		}
		
		++i_pos;
		++j_pos;
		
	}
	// atualizando numero de parametros
	_NFitPar = npar;
	_NPrivPar = npriv;
	
	// Atualizando parametros para funcao de minimizacao
	
	gsl_vector_free(x);
	gsl_vector_free(ss);
	gsl_multimin_fminimizer_free(s);

	x = gsl_vector_alloc(_NFitPar+_NPrivPar*(nData()-1));
	ss = gsl_vector_alloc(_NFitPar+_NPrivPar*(nData()-1));
	minex_func.n = _NFitPar+_NPrivPar*(nData()-1);
	
	for (int i = 0, j = 0; i < _NFitPar; ++i){
		
		if (_PrivParams[_FitPar_index[0][i]][_FitPar_index[1][i]] > 0.0){
			// Parametro eh privado de cada curva de luz
			for (int k = 0; k < nData(); ++k) {
				gsl_vector_set(x, i+k, getpar(_FitPar_index[0][i], _FitPar_index[1][i]));
				gsl_vector_set(ss,i+k, _StepSize[i]); // tamanho do passo p/ funcao de minimizacao
				
			}
			j += nData()-1;
			
		}
		else {
			gsl_vector_set(x, i+j, getpar(_FitPar_index[0][i], _FitPar_index[1][i])); // Parametro para funcao de minimizacao
			gsl_vector_set(ss, i+j, _StepSize[i]); // tamanho do passo p/ funcao de minimizacao
		}
		
		
	}
	minex_func.n = _NFitPar+_NPrivPar*(nData()-1);
	minex_func.f = &chi_square;
	minex_func.params = (void *)this;
	
	init();
}

void inline GSL_Simplex::init(void){
	s = gsl_multimin_fminimizer_alloc (T, _NFitPar+_NPrivPar*(nData()-1));
	gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

}

// Fim da fucao que muda parametros a serem ajustados
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Define como lidar com os parametros em casos extremos!

double inline GSL_Simplex::dealer(const int index,double val)
{
	
	if ( ( _FitPar_index[0][index] == 1 ) && ( _FitPar_index[1][index] == 0 ) ) {
		// inclinacao
		return asin(sin( fabs(val) * M_PI / 180.0 )) * 180.0 / M_PI;
	}
	if ( ( _FitPar_index[0][index] == 1 ) && ( _FitPar_index[1][index] == 1 ) ) {
		// M1 > 0
		return val > 0.0 ? val : -val;
	}	
	if ( ( _FitPar_index[0][index] == 1 ) && ( _FitPar_index[1][index] == 2 ) ) {
		// 0 < M2 < M1
		double M1 = getpar(1, 1);
		//double res = val+step > 0.0 ? (val+step) : (val-step);
		//res = res < M1 ? res : 2.0 * M1 - res; 
		double res = val;
		return res < M1 ? res : M1 - fabs( res/M1 - ceil(res/M1) );
		
	}		
	if ( ( _FitPar_index[0][index] == 2 || _FitPar_index[0][index] == 3 ) && ( _FitPar_index[1][index] == 0 ) ) {
		// flux1 e flux2
		return fabs(val); // > 0.0 ? (val+step) : (val-step);
	}		
	if ( ( _FitPar_index[0][index] == 2 || _FitPar_index[0][index] == 3 ) && ( _FitPar_index[1][index] == 5 ) ) {
		// ffac1 e ffac2
		double rval = val >= 1.0 ? (val) : (1.0+fabs(val));
		
		return rval < 200. ? fabs(rval) : 200.0 - fabs(val/200. - ceil(val/200.));
	}
	else {
		return val;
	}
	
}

//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Calcula o X2 dos dados em relacao ao modelo

double inline GSL_Simplex::Evaluate(const gsl_vector *par)
{

	double* X2 = new double [nData()];
	double* X2_iter = X2;

	for (int ncl = 0; ncl < nData(); ++ncl) {
	
		// Colocar parametros onde deveriam 
		
		for (int np = 0; np < minex_func.n ; ++np) {
			
			
			if (  ( (_FitPar_index[0][np] == _PrivPar_index[0][np]) && (_FitPar_index[1][np] == _PrivPar_index[1][np]) ) ) {
			
				setpar(_FitPar_index[0][np], _FitPar_index[1][np], dealer(np,gsl_vector_get(par, np+ncl)));
				//par+=nData();
			
			}
			else {
				setpar(_FitPar_index[0][np], _FitPar_index[1][np], dealer(np,gsl_vector_get(par, np)));
			}

			
		}
		
		Roche star_2(slaces(), org_par(0));
		Roche star_1(slaces(), org_par(1));
	
		*X2_iter = 0.0;
		
		for (int nphase = 0; nphase < nData(ncl); ++nphase) {
			
			double fmod = star_2.flux<Roche>(DataReaderFase(ncl,nphase),star_1,oculta) + star_1.flux<Roche>(DataReaderFase(ncl,nphase)-0.5,star_2,oculta);
			
			*X2_iter += (DataReaderFlux(ncl,nphase) - fmod) * (DataReaderFlux(ncl,nphase) - fmod) / DataReaderFErr(ncl,nphase);
			
			//cout << DataReaderFase(ncl,nphase) << " " << DataReaderFlux(ncl,nphase) << " " << fmod << " " << DataReaderFErr(ncl,nphase) << endl;
			
		}

		//*X2_iter /= (nData(ncl) - minex_func.n);
		++X2_iter;
		
	}
	
	double res = 0.0;
	X2_iter = X2;
	for (int i = 0; i < nData(); ++i, ++X2_iter) res+=*X2_iter;
	
	res /= double(nData());
	
	delete [] X2;
	
	return res;
	 
}


//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// operador pre incremento (++Object) define iteracao

int GSL_Simplex::operator++ ()
{
	++_iter;
	
	_status = gsl_multimin_fminimizer_iterate(s);
	
	return _iter;
}
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Guarda resultados finais

void GSL_Simplex::store_res(void)
{
	ostringstream saida0;
	
	saida0 << _OutputRootFilename << "_" << _iter << ".dat";
	
	ofstream oo(saida0.str().c_str());
		
	for (int i = 0; i < minex_func.n; ++i) {
		oo << gsl_vector_get(s->x, i) << " ";
	}
	oo << endl;

	oo.close();
	
	cout << "# - Resultdado salvo em: " << saida0.str() << endl;
	
}

//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Guarda curva de luz

void GSL_Simplex::store_cl(void)
{
	for (int ncl = 0; ncl < nData(); ++ncl) {
		
		// Colocar parametros onde deveriam 
		
		for (int np = 0; np < minex_func.n ; ++np) {
			
			
			if (  ( (_FitPar_index[0][np] == _PrivPar_index[0][np]) && (_FitPar_index[1][np] == _PrivPar_index[1][np]) ) ) {
				
				setpar(_FitPar_index[0][np], _FitPar_index[1][np], dealer(np,gsl_vector_get(s->x+ncl, np)));
				//par+=nData();
				
			}
			else {
				setpar(_FitPar_index[0][np], _FitPar_index[1][np], dealer(np,gsl_vector_get(s->x, np)));
			}
			
			
		}
		
		Roche star_2(slaces(), org_par(0));
		Roche star_1(slaces(), org_par(1));
		
		ostringstream saida0;
		
		saida0 << _OutputRootFilename << "_cl" << ncl << "_" << _iter << ".dat";
		
		ofstream oo(saida0.str().c_str());
		
				
		for (int nphase = 0; nphase < nData(ncl); ++nphase) {
			
			double fmod = star_2.flux<Roche>(DataReaderFase(ncl,nphase),star_1,oculta) + star_1.flux<Roche>(DataReaderFase(ncl,nphase)-0.5,star_2,oculta);
			
			//*X2_iter += (DataReaderFlux(ncl,nphase) - fmod) * (DataReaderFlux(ncl,nphase) - fmod) / DataReaderFErr(ncl,nphase);
			
			oo << DataReaderFase(ncl,nphase) << " " << fmod << " " << DataReaderFlux(ncl,nphase) << " " << DataReaderFErr(ncl,nphase) << endl;
			
		}
		oo.close();
		cout << "# - Curva de luz " << ncl << " salvo em: " << saida0.str() <<endl;
	}
	
}


//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
//
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Define funcao chi^2

double chi_square(const gsl_vector *par, void * GSLSimplex)
{
	GSL_Simplex* Simplex_run = (GSL_Simplex *) GSLSimplex;

	return Simplex_run->Evaluate(par);
}
